args <- commandArgs(TRUE)
if(length(args) != 3){
	cat ("Usage:\n\tRscript Q20Q30.R <Distribution_of_Q20_Q30_bases_by_read_position_1.txt> <Distribution_of_Q20_Q30_bases_by_read_position_2.txt> <output png file>\n")
	q(save="no")
}
rt1_1 <- read.table(pipe(paste("cut -f1-3 ",args[[1]])),skip=1)
rt1_2 <- read.table(pipe(paste("cut -f4,5 ",args[[1]])),skip=1)
rt2_1 <- read.table(pipe(paste("cut -f1-3 ",args[[2]])),skip=1)
rt2_2 <- read.table(pipe(paste("cut -f4,5 ",args[[2]])),skip=1)
pos <- c(rt1_1$V1,rt2_1$V1+max(rt1_1$V1))
rawq20 <- c(as.numeric(sub("%","",rt1_1$V2)),as.numeric(sub("%","",rt2_1$V2)))
rawq30 <- c(as.numeric(sub("%","",rt1_1$V3)),as.numeric(sub("%","",rt2_1$V3)))
cleanq20_1 <- as.numeric(sub("%","",rt1_2$V1))
cleanq20_2 <- as.numeric(sub("%","",rt2_2$V1))
cleanq30_1 <- as.numeric(sub("%","",rt1_2$V2))
cleanq30_2 <- as.numeric(sub("%","",rt2_2$V2))
png(file=args[[3]],height=360,width=576,antialias='default')
plot(pos,rawq20,main="Q20 Q30 base percentage along reads",type="l",col="red",lwd=2, bty="o",xaxt="n",yaxt="n",xlab="",ylab="",xaxs="i",ylim=c(0,100))
lines(pos,rawq30,lty=1,col="green",lwd=2)
pos1 <- seq(length(rt1_1$V1)-length(rt1_2$V1)+1,length(rt1_1$V1))
pos2 <- seq(length(rt2_1$V1)-length(rt2_2$V1)+1,length(rt2_1$V1))+max(rt1_1$V1)
lines(pos1,cleanq20_1,lty=2,col="red",lwd=2)
lines(pos2,cleanq20_2,lty=2,col="red",lwd=2)
lines(pos1,cleanq30_1,lty=2,col="green",lwd=2)
lines(pos2,cleanq30_2,lty=2,col="green",lwd=2)
lines(c(max(rt1_1$V1),max(rt1_1$V1)),c(0,100),col="blue",lty=4,lwd=2)
axis(side=1, seq(0,max(pos),20), tcl=0.2, labels=FALSE)
axis(side=2, seq(0,100,20), tcl=0.2, labels=FALSE)
axis(side=3, seq(0,max(pos),20), tcl=0.2, labels=FALSE)
axis(side=4, seq(0,100,20), tcl=0.2, labels=FALSE)
mtext(seq(0,max(pos),20),side=1,las=1,at=seq(0,max(pos),20),line=0.3,cex=1.5)
mtext(seq(0,100,20),side=2,las=1,at=seq(0,100,20),line=0.3,cex=1.5)
mtext("Position along reads",side=1, line=2, at=max(pos)/2, cex=1.5)
mtext("Percent",side=2, line=2.5, at=50, cex=1.5)
legend("bottomleft",lwd=2,lty=c(1,1,2,2),col=c("red","green","red","green"),bty="n",legend=c("raw Q20","raw Q30","clean Q20","clean Q30"),cex=1,inset=0.05)
d <- dev.off()
